#include "../mds/frame.h"

static 
void mds_frame_reset( mds_frame *frame )
{
	frame->num_particles = 0;
	frame->num_sets      = 0;
	frame->sets          = NULL;
	frame->reserved      = 0;
	frame->nodes         = NULL;
}

int  mds_frame_load( mds_frame *frame, size_t num_particles, size_t num_sets )
{
	assert( frame != NULL );
	assert( num_particles > 0 );
	assert( num_sets > 0 );
	
	mds_frame_reset( frame );
	
	
	frame->sets = MDS_CALLOC_AS(mds_vector_set,num_sets);
	if( NULL == frame->sets )
	{
		return -1;
	}
	
	frame->num_particles = num_particles;
	frame->num_sets      = num_sets;
	
	while( frame->reserved < frame->num_sets )
	{
		if( NULL == (frame->sets[ frame->reserved ] = mds_vector_set_create( frame->num_particles ) ) )
		{
			mds_frame_quit( frame );
			return -1;
		}
		frame->reserved++;
	}
	
	frame->nodes = mds_pnode_create(num_particles);
	if( NULL == frame->nodes )
	{
		mds_frame_quit(frame);
		return -1;
	}
	
	return 0;
	
}

void mds_frame_quit( mds_frame *frame )
{
	assert( frame );
	
	mds_pnode_delete( frame->nodes, frame->num_particles );
	while( frame->reserved > 0 )
	{
		mds_vector_set_delete( frame->sets[ --(frame->reserved) ], frame->num_particles);
	}
	MDS_FREE_AS(mds_vector_set, frame->sets, frame->num_sets );
	
	mds_frame_reset(frame);
}


void mds_frame_initialize_positions( mds_frame *frame )
{
	assert( frame );
	assert( frame->pbc_for_r );
	assert( frame->sets      );
	assert( frame->num_sets > MDS_POSITION );
	mds_position_initialize( frame->sets[MDS_POSITION], frame->num_particles, frame->pbc_for_r, & frame->box );
}

void mds_frame_initialize_velocities( mds_frame *frame, Real Temperature)
{
	assert( frame );
	assert( frame->sets      );
	assert( frame->num_sets >= MDS_VELOCITY );
	frame->K = mds_velocity_initialize( frame->sets[MDS_POSITION], frame->sets[MDS_VELOCITY], frame->num_particles, Temperature );
	
}

void mds_frame_update_kinetic_energy( mds_frame *frame )
{
	assert( frame );
	assert( frame->sets      );
	assert( frame->num_sets >= MDS_VELOCITY );
	
	frame->K = mds_velocity_kinetic_energy( frame->sets[MDS_VELOCITY], frame->num_particles);
}

void mds_frame_no_forces( mds_frame *frame )
{
	assert( frame );
	assert( frame->sets      );
	assert( frame->num_sets >= MDS_ACCELERATION );
	{
		register size_t i;
		mds_vector     *f = frame->sets[ MDS_ACCELERATION ];
		for( i= frame->num_particles; i >0 ; --i, ++f )
		{
			mds_vector_ldz( f );
		}
	}
	
}

void mds_frame_rescale( mds_frame *frame, const mds_vector *box )
{
	mds_vector  scal = MDS_VECTOR_INIT;
	assert( frame );
	assert( frame->num_particles > 0 ); 
	assert( box != NULL );
	assert( frame->pbc_for_r != NULL );
	assert( frame->num_sets > 0 );
	
	scal.x = box->x / frame->box.x;
	scal.y = box->y / frame->box.y;
	scal.z = box->z / frame->box.z;
	fprintf( stderr, "@@ rescaling by [%g %g %g]\n", scal.x, scal.y, scal.z);
	
	/** TODO: change according to PBC ? **/
	
	{
		const size_t  np  = frame->num_particles;
		mds_vector   *pos = frame->sets[MDS_POSITION];
		size_t       i;
		for( i=0; i < np; ++i )
		{
			mds_vector *r = &pos[i];
			frame->pbc_for_r(r,&frame->box);
			r->x *= scal.x;
			r->y *= scal.y;
			r->z *= scal.z;
			/* frame->pbc_for_r(r,box); */
		}
		
		mds_vector_cpy( &frame->box, box);
		
	}
	
	
}


#include "./constants.h"
Real mds_frame_T( const mds_frame *frame )
{
	return frame->K / ( MDS_F(0.5) * ( MDS_F(3.0) * frame->num_particles) * MDS_R );
}

void mds_frame_NVE_T_and_PV( const mds_frame *frame, Real *T, Real *PV)
{
	assert( frame );
	assert( T );
	assert( PV );
	{
		*T  = mds_frame_T(frame);
		*PV = frame->num_particles * MDS_R * (*T) + frame->W; 
	}
	
}

